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Abstract 

The parallel computational complexity of the Bak-Sneppen evolution model is studied. It is 
shown that Bak-Sneppen histories can be generated by a massively parallel computer in a time 
that is polylogarithmic in the length of the history. In this parallel dynamics, histories are built 
up via a nested hierarchy of avalanches. Stated in another way, the main result is that the logical 
depth of producing a Bak-Sneppen history is exponentially less than the length of the history. This 
finding is surprising because the self-organized critical state of the Bak-Sneppen model has long 
range correlations in time and space that appear to imply that the dynamics is sequential and 
history dependent. The parallel dynamics for generating Bak-Sneppen histories is contrasted to 
standard Bak-Sneppen dynamics. Standard dynamics and an alternate method for generating his- 
tories, conditional dynamics, are both shown to be related to P-complete natural decision problems 
implying that they cannot be efficiently implemented in parallel. 
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I. INTRODUCTION 



The Bak-Sneppen model [|], Q is among the simplest systems to display self-organized crit- 
icality and avalanche behavior. Starting from arbitrary initial conditions, the Bak-Sneppen 
model spontaneously reaches a state with long range correlations in space and time that are 
set up by a series of avalanche events. The approach to criticality, the critical state of the 
model, and the properties of avalanches in the Bak-Sneppen model have been extensively 
studied during the past decade. There are now precise values of critical exponents and well- 
tested scaling laws. On the other hand, there is not yet a theory of the critical state except 
for the case of mean field theory and near zero dimension. 

In this paper we investigate the complexity of the histories generated by the Bak-Sneppen 
model. These histories manifest long time correlations with early random events influencing 
distant, later events. However, the past may determine the future in ways which are more 
or less complex. It is the purpose of this paper to give a precise meaning to the notion of 
a complex history and then to show that histories produced by the Bak-Sneppen model are 
not complex in this sense. 

The framework in which we will define the complexity of a history is computational 
complexity theory. Computational complexity theory measures the computational resources 
needed to solve problems. Among those resources are parallel time or equivalently depth, 
which measures the number of steps needed by a massively parallel computer to solve a 
problem. If the solution to a problem can be re-organized so that many computations are 
carried out simultaneously and independently then a great deal of computational work can 
be accomplished in few parallel steps. On the other hand, it may be that the sequential 
nature of the problem precludes reorganizing the work into a small number of parallel steps. 
The theory of parallel computational complexity provides algorithms for efficiently solving 
some problems in parallel and provides evidence supporting the proposition that solutions 
to other tasks cannot be substantially accelerated by parallelization. A powerful feature 
of this theory is that it yields essentially the same results independent of the choice of 
the model of computation. Thus, measures of parallel computational complexity reveal the 
logical structure of the problem at hand rather than the limitations of the particular model 
of parallel computation. 

We will define the complexity of histories generated by a stochastic model such as the 
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Bak-Sneppen model in terms of the parallel time required to generate a sample history. 
If, through parallelization, the dynamics of the model can be reorganized so that many 
processors running for a few steps can generate a statistically valid history then we will 
say the histories generated by the model are not complex. If on the other hand, there 
is little benefit from parallelization so that a history must be generated in a step by step 
fashion, then the model generates complex histories. The Bak-Sneppen model illustrates the 
point that scale invariant spatial and temporal correlations are not sufficient to guarantee 
complexity in the sense defined here. Indeed, it is the strictly hierarchical nature of the 
avalanches in the Bak-Sneppen model that permits us to generate them in a small number 
of highly parallel steps. 

The present work fits into a series of studies of the computational complexity of models in 
statistical physics. Although previous studies focused mainly on the computational complex- 
ity of generating spatial patterns, in many cases the analysis can also be applied to reveal the 
complexity of the history of the growth of the pattern. For example, two dynamical processes 
for diffusion limited aggregation (DLA) have been shown to be P-complete 0, These 
results suggest that these processes for generating DLA clusters cannot be re-organized to 
greatly reduce the parallel time needed for the generation of an aggregate (however, see Q 
for a parallel algorithm that achieves a modest speed-up for DLA). The growth of an aggre- 
gate can be viewed as a history and the same arguments can be employed to give evidence 
that such histories cannot be significantly parallelized. Similarly, predicting the behavior of 
Abelian sandpiles for d > 3 has been shown to be P-complete By contrast, the patterns 
generated by a number of growth models such as the Eden model, invasion percolation and 
the restricted solid-on-solid model were shown to require little parallel time [0 and these 
results carry over immediately to fast parallel algorithms for generating histories for these 
models. 

Summary 

In Sec. II we give an overview of selected aspects of computational complexity theory 
relevant to the present work and in Sec. Ill we review the Bak-Sneppen model. In Sec. 
we present parallel dynamics for the Bak-Sneppen model and in Sec. we present simula- 
tion results showing that Bak-Sneppen histories can be efficiently simulated using parallel 
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dynamics. In Sec. ^ we show that the standard dynamics for the Bak-Sneppen model and 
a new dynamics based on conditional probabilities both yield P-complete problems. The 
paper concludes in Sec. |V11. 



II. COMPUTATIONAL COMPLEXITY 

Computational complexity theory determines the scaling of computational resources 
needed to solve problems as a function of the size of the problem. Introductions to the 
field can be found in Refs. . Computational complexity theory may be set up in sev- 

eral nearly equivalent ways depending on the choice of model of computation. Here we focus 
on parallel computation and choose the standard parallel random access machine (PRAM) 
as the model of computation. The main resources of interest are parallel time or depth and 
number of processors. A PRAM consists of a number of simple processors (random access 
machines or RAMs) all connected to a global memory. Although a RAM is typically defined 
with much less computational power than a real microprocessor such as Pentium, it would 
not change the scaling found here to think of a PRAM as being composed of many micropro- 
cessors all connected to the same random access memory. The processors run synchronously 
and each processor runs the same program. Processors have an integer label so that different 
processors follow different computational paths. Since all processor access the same memory, 
provision must be made for how to handle conflicts. Here we choose the concurrent read, 
concurrent write (CRCW) PRAM model in which many processors may attempt to write 
to or read from the same memory cell at the same time. If there are conflicts between what 
is to be written, the processor with the smallest label succeeds in writing. 

The PRAM is an idealized and fully scalable model of computation. The number of 
processors and memory is allowed to increase polynomially (i.e. as an arbitrary power) in 
the size of the problem to be solved. Communication is idealized in that it is assumed that 
any processor can communicate with any memory cell in a single time step. Obviously, this 
assumption runs up against speed of light or hardware density limitations. Nonetheless, 
parallel time on a PRAM quantifies a fundamental aspect of computation. During each 
parallel step, processor must carry out independent actions. Communication between pro- 
cessor can only occur from one step to the next via reading and writing to memory. Any 
problem that can be solved by a PRAM with H processors in parallel time T could also be 
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solved by a single processor machine in a time that is no greater than HT since the single 
processor could sequentially run through the tasks that were originally assigned to the H 
processors. On the other hand, it is not obvious whether the work of a single processor can 
be re-organized so that it can be accomplished in a substantially smaller number of steps by 
many processor working independently during each step. Two examples will help illustrate 
this point. First, suppose the task is to add n finite-precision numbers. This can be done 
by a single processor in a time that scales linearly in n. On a PRAM with n/2 processor 
this can be done in O(logn) parallel time using a binary tree (for simplicity, suppose n is a 
power of 2). In the first step, processor one adds the first and second numbers and puts the 
result in memory, processor two adds the third and fourth numbers and puts the result in 
memory and so on. After the first step is concluded there are n/2 numbers to add and these 
are again summed in a pairwise fashion by n/4 processors. The summation is completed 
after O(logn) steps. This method is rather general and applies to any associative binary 
operation. Problems of this type have efficient parallel algorithms and can be solved in time 
that is a power of the logarithm of the problem size, here n, that is, polylog time. 

A second example is the problem of evaluating the output of a Boolean circuit with 
definite inputs. A Boolean circuit is composed of AND, OR and NOT gates connected by 
wires. Suppose the gates are arranged in levels so that gates in one level take their inputs 
only from gates of the previous level and give their outputs only to gates of the next level 
so there is no feedback in the circuit. The number of levels in the circuit is referred to as 
the depth of the circuit. At the top level of the circuit are the TRUE or FALSE inputs. 
At the bottom level of the circuit are one or more outputs. Given some concisely encoded 
description of the circuit and its inputs, the problem of obtaining the outputs is known as 
the circuit value problem (CVP). Clearly, we can solve CVP on a PRAM in parallel time that 
is proportional to the depth of the circuit since we can evaluate a single level of the circuit 
in a single parallel time step. On the other, there is no known general way of speeding up 
the evaluation of a Boolean circuit to a time that is, say polylog in the depth of the circuit. 
For the case of adding n numbers, the logical structure of the problem is sufficiently simple 
that we can substitute hardware for time to obtain the answer in polylog parallel steps. For 
CVP, hardware cannot be used to substantially decrease the depth of the problem. 

At the present time, the statement that CVP cannot be solved in polylog parallel time is 
not mathematically proved. Instead, one can show that CVP is P-complete. To understand 
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the meaning of P-completeness we must first introduce the complexity classes P and NC 
and the notion of reduction. P consists of the class of problems that can be solved by a 
PRAM in polynomial time with polynomially many processors. NC consists of the class of 
problems that can be solved in polylog time on a PRAM with polynomially many processors. 
A problem A is reduced to a problem i3 if a PRAM with an oracle for B can be used to solve 
A in polylog time with polynomially many processors. An oracle for B is able to supply 
the PRAM a solution to an instance of i3 in a single time step. Intuitively, if A can be 
reduced to B then A is no harder to solve in parallel than B. A problem in P is P-complete 
if all other problems in P can be reduced to it. It can be proved that CVP is P-complete. 
Furthermore, if it could be shown that there is any problem in P that is not in NC then 
it would follow that CVP is not in NC . However, the conjecture that P 7^ NC, though 
widely believed, is not yet proved. Since reductions are transitive, showing that another 
problem A is P-complete can be proved by reducing CVP to A. Later in the paper we will 
take this approach to prove that two problems representing different Bak-Sneppen dynamics 
are P-complete. 

In addition to supplying a P-complete problem. Boolean circuits serve as an alternate 
model of computation. For a given computational problem, and for a given problem, there 
is a Boolean circuit that solves the problem. In order to solve problems of any size, a family 
of Boolean circuits is needed, one for each problem size. Furthermore, to be comparable 
to a PRAM running a single program, the family of circuits should satisfy a uniformity 
condition that bounds the complexity of the computation needed to specify each circuit in 
the family. The resources of circuit depth (number of levels in the circuit) and circuit width 
(maximum number of gates in a level) are roughly comparable to the PRAM resources of 
parallel time and number of processor, respectively. Indeed the class NC can be defined 
as the set of problems that can be solved by uniform circuit families whose depth scales 
polylogarithmically in the problem size and whose width scales polynomially in the problem 
size. P is the class of problems solvable by circuits of polynomial width and depth. 

III. BAK-SNEPPEN MODEL 

In this section we review the Bak-Sneppen model and its behavior. A detailed discussion 
can be found in Ref. The Bak-Sneppen model defines histories of configurations on a d- 
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dimensional lattice. At each lattice site i and discrete time t, there is a number Xi{t) G [0, 1]. 
The standard dynamics for the Bak-Sneppen model creates a history sequentially with one 
new time slice determined in each step of the dynamics. A single time step in standard 
dynamics for the Bak-Sneppen model consists of locating the extremal site s, which has 
the minimum number. The extremal site and its 2d nearest neighbors are given new x's 
chosen randomly from the uniform distribution on [0, 1]. The x's on the remaining sites are 
unchanged. Let f{t) be the extremal value at the extremal site s{t) at time t. 

Bak-Sneppen histories can be organized into avalanches or causally connected regions 
of activity in the d + 1-dimensional space-time. Avalanches can be classified by bounding 
the extremal value during the avalanche. The collection of renewed sites from t to t' is an 
/-avalanche if f{t - 1) > / and f{t') > f but for all t" with t < t" < t', f{t") < f. Note 
that it is also possible to have an /-avalanche whose duration is a single time step (i. e. 
t = t' and no times t") if f{t — 1) > / and f{t) > f. An /-avalanche can be partitioned into 
a sequence of one or more /'-avalanches for f < f. Repeated use of this idea allows us to 
partition avalanches into a nested sequence of subavalanches. The partitioning of avalanches 
into subavalanches is one of the key ideas in the parallel algorithm described below. 

The properties of avalanches are entirely independent of the environment in which they 
exist so long as the environment is large enough to contain the avalanche. Indeed we can 
grow the backbone of an avalanche without regard to the environment. The backbone of an 
/-avalanche is the set of sites i and the corresponding Xi for which Xi < f. (Note that the 
last time step in an avalanche is not part of the backbone). The backbone of an /-avalanche 
for / = 0.5 is shown in Fig. ^ The only information we need about the environment is that 
all environmental sites have random numbers uniformly distributed in the interval [/, 1]. 
The values of these numbers is unimportant and the full lattice of random numbers does 
not have to be realized. During the construction of an /-avalanche, sites must be renewed. 
Instead of assigning random numbers in the interval [0, 1] as is done in standard Bak-Sneppen 
dynamics, one can declare a renewed site to be part of the environment with probability 
1 — / or, with probability /, declare it part of the backbone. If the site is part of the 
backbone it is then given an explicit random number in the interval [0, /). The construction 
of an /-avalanche backbone begins with the renewal of the origin and its neighbors and then 
continues until there are no remaining sites with number less than /. This approach for 
making an avalanche backbone is called the Bak-Sneppen branching process in |]T0| and is 
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an efficient method for carrying out large-scale simulations of the Bak-Sneppen model. 
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FIG. 1: The backbone of an /-avalanche for / = 0.5 in the one-dimensional Bak-Sneppen model. 
Black dots represent positions where the random numbers are less than /. Immediately before and 
at the last time of the avalanche all random numbers are greater than /. 

The Bak-Sneppen model manifests self-organized criticality: independent of the initial 
state, after a long time the system approaches a stationary state where the average avalanche 
size diverges and the distribution of avalanche sizes follows a power law. In the critical state 
almost all x's exceed a critical value, fc (for the = 1, ~ 0.6670). If S is the duration of 
an avalanche then, in the critical state, the probability of an avalanche of duration S decays 
as S~'^ (r ^ 1.1 for d = 1). The average duration {S)j of /-avalanches diverges as 

{s)f - {fc - fY' (1) 

with 7 !^ 2.7 for d = 1. The /-avalanche size distribution, V{S, /) has the scaling form 
V{S, f) = S^'^g{S{fc — fY^'^) with a, r and 7 related by the usual exponent relation 
7 = (2 — r)/cr. The scaling function g{y) goes to a constant as y ^ and decays rapidly 
to zero for y ^ 1. Thus the upper cut-off in the /-avalanche size distribution obeys the 
scaling law Sco ~ {fc — fY^"^ (o" ~ 0.34 for ci = 1 |[ll|]). The spatial extent of /-avalanches 
also diverges in the critical state. Let W be the number of sites covered by an /-avalanche, 
then {W)f = {f,~f)-\ 

The approach to the critical state is characterized by a gap G{t), which is a stepwise 
increasing function of time and is the maximum value of f{t') for t' < t. Times when the 
gap increases mark the end of an avalanche. If the gap first has the value G at time t then 
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all x's in the system are greater than or equal to G at time t and a G-avalanche begins at 
time t + 1. On average, the gap approaches the critical value as a power law, 

/,-(G(t))~(t/iV)^ (2) 

where N is the number of sites in the system. In the parallel dynamics for the Bak-Sneppen 
model we will first identify the successive values of the gap and then fill in the intervals 
between them with the appropriate avalanches. 



IV. PARALLEL BAK-SNEPPEN DYNAMICS 



In this section we present a parallel approach for constructing Bak-Sneppen histories. 
At the core of this approach is a procedure for hierarchically building avalanches in which 
small avalanches are pasted together to yield larger avalanches. The large scale structure of 
the history is controlled by the sequence of gap values. Each value of the gap determines 
an avalanche that needs to be constructed. The algorithm proceeds in three stages. In the 
first stage, described in Sec. [IV A| , the gap sequence is obtained and in the second stage, 



described in Sec. [IV B|, the required avalanches are constructed. Finally, in the third stage, 
all space-time points that are not yet determined are explicitly given numbers. 

Although we have not run the parallel algorithm on a parallel computer, we have measured 
essential aspects of its behavior using a conventional Bak-Sneppen simulation and the results 
are reported in Sec. 0. The key finding is that we can produce Bak-Sneppen histories in a 
parallel time that is polylogarithmic in the length of the history. 



A. Parallel algorithm for the gap sequence 

In this subsection we show how to construct a sequence of gap values, Gk-, k = 1, . . . , K 
for a system with lattice sites (independent of dimension). Suppose that at some time 
the gap increases from Gk to Gfc+i- This event occurs at the last step of a G/j-avalanche so, 
before the selection of the extremal value that increases the gap, the distribution on every 
site of the lattice is uniform on the interval [Gk, !]• Thus Gk+i is a random variable given 
by 

Gfc+i = Gk + X{1 -Gk) = Gk(BX (3) 
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where X is a random variable that is the minimum over N random numbers uniformly 
distributed on [0, 1] and © is the XOR operation for real numbers. Since © is associative 
we can unroll Eq. (||) to obtain, 

Gk = Go®Xi®X2®---®Xk (4) 

where Go is the initial value of the gap and the Xj are identically distributed independent 
random variables chosen from the same distribution as X described above. Using a standard 
parallel minimum program, each Xj can be generated on a probabilistic PRAM in 0{\ogN) 
parallel time. A standard parallel prefix algorithm [^] can then be used to evaluate all K 
values of the gap in 0{logK + logN) parallel time. 

B. Parallel algorithm for /-avalanches 

In this section we show how to generate an /-avalanche in parallel time that is polylog in 
ifc — /)• The idea of this algorithm is to generate in parallel a large number of independent 
small avalanches and then to paste groups of these together to create a smaller number of 
larger avalanches. These are in turn pasted together to create yet larger avalanches and so 
on in a hierarchical fashion. The algorithm is controlled by a pre-determined list of / values, 
/o < /i < • • • < /• These values specify the types of avalanches to be manufactured. The 
sequence of / converges to fc as 

/„ = /,(1 - (5) 

with q a constant between to 1. This choice for /„ implies that for large n, {S) j^/ {S) f^_^ —>■ 
q~"' so that the average number of /„_i-avalanches needed to construct an /„-avalanche is 
asymptotically constant (and equal to q~^). The length of the sequence and its penultimate 
element /„^^^ is chosen so that /„„,^^ < / < fn^^^+i- 

Since /o = 0, /o-avalanches always have S = 1 and consist of the origin and its 2d nearest 
neighbors together with their associated x values. In parallel we can make M, /o-avalanches 
in constant time using M processors. We now construct /i-avalanches by pasting together 
one or more /o-avalanches. The origin of the first /o-avalanche is also the origin of the 
/i-avalanche. If the minimum of the 2d + 1 x's of the /o-avalanche is greater than /i then 
we are done and the /o-avalanche is also an /i-avalanche. If not, a second /o-avalanche is 
attached to the first /o-avalanche with its origin lined up with the extremal site of the first 
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avalanche. Again, we check whether the minimum x among all the sites covered by the pair 
of /o-avalanches is greater than /i and, if so, we are done. The pattern for constructing /„- 
avalanches from an /„_i-avalanche should now be clear: /„_i-avalanches are pasted together 
until all the covered sites have a; > /„. 

The process of pasting together /„_i-avalanches is a generalized version of the standard 
Bak-Sneppen process. Indeed, pasting together /o-avalanches is exactly the standard Bak- 
Sneppen process. All that we need to know about each /„_i-avalanche is its origin and the 
final values of x on the sites the avalanche covers. During the construction of an /„-avalanche, 
when a new /n_i-avalanche is attached, its origin is lined up with the extremal site among all 
of the sites currently covered by the growing /„-avalanche. The x vahics on the covered sites 
of the newly added /„_i-avalanchc renew the sites covered by this avalanche and, if some of 
these sites were previously part of the environment, the total number of covered sites of the 
growing /^-avalanche is increased. If the covered sites all have x > fn, the /„-avalanche is 
finished. The pasting together of /„_i-avalanches must be carried out sequentially so that if 
the number of /n_i-avalanches needed to construct a single /^-avalanches is Kn then it will 
take Kn parallel steps to grow the avalanche. Finally, once an /-avalanche is constructed, it 
is straightforward to find its backbone by eliminating those sites for which Xi > f. 

On average g"''' /„_i-avalanches must be pasted together to form one /„-avalanche, that 
is {Kn) — q"'^- However, the n*^ step in the parallel construction of an /-avalanche is not 
complete until all the required /^-avalanches have been made. Thus, the time required for 
this step will be set by the maximum value of Kn- As / approaches /c, and for fixed n 
we need to make increasingly many /„-avalanches so that the tail of the Kn distribution 
is explored. As we shall see in the next section, the expected value of the maximum of 
increases logarithmically in the number of /„+i-avalanches that are constructed. As a result, 
the average parallel time required to generate an /-avalanche is predicted to be proportional 

to l0g^(/c-/). 

The discussion thus far has posited an effectively infinite environment for the avalanche 
to grow in, however this is not required. The same hierarchical construction can just as well 
be used to make avalanches in periodic boundary conditions in a system of size L. Addi- 
tional complications arise for other boundary conditions due to the breaking of translational 
symmetry and this situation will not be discussed here. 

Finally, we are ready to describe how to construct a full Bak-Sneppen history in a system 
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with sites and periodic boundary conditions starting from the initial condition, Go = 
0. First, following the procedure of Sec. [IV A| we construct a sequence of K gap values, 
Go) Gil • • • ) Gk- Then, using the parallel algorithm described above, for each gap value, G^, 
we construct the backbone of a G^-avalanche for all k < K. The site i where Xi = G^ is 
chosen randomly from among all the sites of the lattice and this site serves as the spatial 
origin for the G^-avalanche. The starting time + 1 of the G^-avalanche is obtained by 
summing over the durations of the preceding avalanches, tk = Y.t=o where Si is the 
duration of the Gravalanche. Successive avalanche backbones are concatenated to form 
the backbone of the history which includes all the avalanche backbones and the intervening 
extremal sites that have the gap values. 

Thus far only the backbone of the history is specified while sites in the environment are 
conditioned by the extremal values but not fully specified. A definite, statistically valid, 
history is constructed from the backbone by fixing all space-time points Xi{t) not in the 
backbone. At the last time, tmax, the gap has the values Gk- The randomly chosen extremal 
site, s(tmax) is part of the backbone and takes the value Gk- All other sites at time tmax are 
chosen from the uniform distribution on [Gk, !]• From these values we now work backwards 
in time independently for each site. Consider site i, we have that Xi{t) = Xj(tinax) until 
the most recent time that i was part of the backbone. During the time intervals that site 
i is part of the backbone, its value is determined. For each connected time interval, [t',t"] 
that site i is part of the environment, Xi takes a single value. There are two possibilities. 
The interval may come to end because i is an extremal site where the gap increases from 
Gk-i to Gk that is, t" = tk and i = s{tk). In this = Gk throughout the interval. 

The second possibility is that site i is renewed at time t" + 1 and becomes part of the G^- 
avalanche backbone at that time. In this chosen uniformly from the distribution 

[Gk, 1] and has this value during the interval from t' to t". Working backwards in this way a 
definite Bak-Sneppen history is re-constructed from the backbone. It should be noted that 
the construction of a definite history from the backbone can be carried out in parallel in 
polylog time since each site can be treated independently and since intervals during which 
a site is not in the backbone can be identified using a parallel graph connectivity algorithm. 
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V. SIMULATION OF PARALLEL DYNAMICS 



We have studied the performance of the parallel algorithm for constructing avalanches in 
the one-dimensional Bak-Sneppen model. We collected statistics on the number of 
avalanches needed in the construction of /„- avalanches during the manufacture of a single 
/„^^^- avalanche. Let fln{n^na.x) = max{K„} be the maximum value of tin required during the 
construction of the /n^^^-avalanche. The parallel running time T(nmax) , up to a constant 
of proportionality, for constructing an /^^^^^-avalanche is obtained by summing On(nniax), 

^^max 

T(nmax) = X! ^n(^max)- (6) 
n=l 

The simulation works by running standard Bak-Sneppen dynamics and, on the fly, keeping 
track of Kn and fin- 

The following pseudo-code sketches a single step of an algorithm that collects statistics 
on the distribution of «;„ and computes the running time for constructing /^^^g^-avalanches. 
It is based on the observation that an /„-avalanche is concluded when the extremal value 
fit) > fn- Arrays «;„ and Q„ are initialized to zero. 

1. Execute one time step of standard Bak-Sneppen dynamics and find the extremal value, 

fit)- 

2. For n — lton — ^max — 1 do 
if f{t) > fn then 

• if > fin then fin <— Kn 

• Collect statistics for 

• Kn ^0 

end do 

3. if f{t) > /nniax then the /„^3^-avalanche ends. Calculate the parallel time, T by 
summing fin- 

else go to Step 1. 
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The re-organization of the sequential history of an avalanche into a hierarchy of sub- 
avalanches is illustrated in Figs. ^ and ^ for the case of an /y-avalanche. In Fig. ^ the sequence 
of extremal values is shown together with horizontal lines indicating bounds fi, . . . , fj. Fig- 
ure ^ shows a tree representing the hierarchy of subavalanches making up the /y avalanche. 
The fn level of the tree represents the /^-avalanches needed in the construction of the final 
/7-avalanche. The number next to each node represents the ending time of a subavalanche. 
For example, the /s level of the tree shows that the /7-avalanche is composed of three /s- 
avalanches ending at times 3, 8 and 9. The total parallel time is the sum of the maximum 
degree of the tree at each level and is 13 in this case. 
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FIG. 2: Extremal values, f{t) vs. time t during an /7-avalanche. Horizontal lines are the values 
/i) • • • ) /?■ 

In our simulations, for each nmax up to 40 we constructed M = 10^ /„^^^-avalanches. The 
sequence of avalanche sizes is given by Eq. |^ with q = 0.9 and fc = 0.667. Figure |^ shows 
the (T(nmax)) vs. n^^^. A fit of the data to the functional form (T(nmax)) = cl^^^x + ^ 
over the range 13 < "n-max < 40 yields z = 1.832 ± 0.004. The fit is reasonable with 
X^/d.o.f = 1.1. We also collected data for q = 0.8 and q = 0.7 over a smaller range of 
n^ax- For q = 0.8 we found that z = 2.00 ± 0.03 (xV^-O-f = 1.3) for 12 < n^ax < 20 
while for q = 0.7 we found that z = 1.88 ± 0.09 (xVd-o.f = 1.3) for 9 < n^^ax < 13. Since 
log(/c — fn) ~ nlogg, a conservative conclusion from the data is that the average parallel 
time (T) for constructing an /-avalanche behaves as (T) = log'^*-^''(/c — /) or, in terms of 
the duration of the avalanche(T) = log^*-^-* S. 
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FIG. 3: Tree representing the re-organization of the /7-avalanche of Fig. |2|into subavalanches. The 
/„-level of the tree represents partitioning the /7-avalanche into /^-avalanches and the numbers on 
the nodes represent the ending times of the /^-avalanches. Roughly speaking, sequential time runs 
to the right and parallel time runs upward. 
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FIG. 4: The average parallel time (T(nmax)) for constructing an /^j^^^^- avalanche vs. ?i-max- The 
straight line is a linear fit to the data for rimax ^ 32. 
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The time requirements of the parallel algorithm can be understood by considering the 
probability density, Pn{K,n) of Recall that «;„ is the number of /„_i-avalanches needed 
in the construction of an /^-avalanche and that its average is = q~'^. Suppose that the 
tail of Pn{K.n) behaves as e"***". Given these assumptions we can estimate the maximum 
value, f^nl^max), of K„ in the construction of an /„_^^^^-avalanche. The estimate is obtained 
by requiring that the expected number of /^-avalanches manufactured times -Pn(^n (^max)) 
is order unity. The expected number of /^-avalanches made during the construction of an 
/rimax'S'Valanche is g-("max-n)7 have 



Summing over n yields 



SJ„(„„„) « - (7) 



r(„_) « !W,M1M (8) 



in qualitative agreement with the numerical results. 

We measured P„(k„) during our simulations. The logarithm of the histogram of k„ for 
several values of n is shown in Fig. |^. It is clear that the behavior of Pn{i^n) is more 
complicated than simple exponential decay. There is a short time fast decay and, for larger 
ra, a long time slow decay. It is possible that the asymptotic behavior of P„(k„) for the two 
largest values of n in Fig. ^both display exponential decay with the same value of h. It is also 
possible that 6 is a slowly decreasing function of n or that the asymptotic decay is slower than 
exponential. If h decreases as a power of n or if the tail of P„ is a stretched exponential it will 
still be the case that (r(nmax)) is polynomial in nmax and thus polylogarithmic in (/c — /). 
On the other hand, if the tail of P„ is a power law then (T(nmax)) would be exponential in 
^max- The data P„ does not suggest a power law tail and the data for (T(?2inax)) appears to 
rule out exponential growth however more numerical work is needed to fully understand the 
asymptotic behavior of these quantities. 



VI. SEQUENTIAL BAK-SNEPPEN DYNAMICS AND P-COMPLETENESS 

We have seen in the previous section that the construction of Bak-Sneppen histories can 
be arranged to be carried out in parallel in a time that is polylog in the length of the 
history. There are however alternative dynamics for producing histories with the correct 
statistical properties. In this section we show that two of these dynamics have the property 
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FIG. 5: Histogram of k„ for several avalanches. 

of P-completenesss. The property of P-completeness imphes that it is not possible to run 
these dynamics in parallel in polylog time. The two dynamics considered here are the 
standard dynamics by which the Bak-Sneppen model is usually defined and a conditional 
dynamics where explicit values are not assigned to sites until it is necessary. Thus, if we 
knew only about standard or conditional dynamics we would conclude that the Bak-Sneppen 
model generates complex histories. Parallel Bak-Sneppen dynamics shares features of both 
standard dynamics and conditional dynamics. 

Although standard Bak-Sneppen dynamics cannot be implemented in parallel in polylog 
time, in Sec. [VI C| we show how to use the parallel construction of avalanches to achieve 
power law speed-up. 

A. Standard Bak-Sneppen dynamics 

Standard Bak-Sneppen dynamics formalizes the usual approach to generating Bak- 
Sneppen histories. Initially, the lattice is described by values {xj(0)}. A list of numbers 
yin) ^ j-Q^ 1], n = 1, 2, ... , {2d + l)tmax} controls the dynamics. At each time step, the 
extremal site is located. The extremal site and its 2d neighbors are renewed using the next 
group oi 2d + 1 y values. Periodic boundary conditions are assumed. Figure ^ shows an 
example oi d = 1 standard Bak-Sneppen dynamics. The numbers along the top row of the 
figure are {xj(0)} and the numbers in the three columns on the right are the groups of y's. 
Successive states of the lattice are given in successive rows in the lower left section of the 
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FIG. 6: An example of standard Bak-Sneppen dynamics. The initial values are above the horizontal 
line and time proceeds downward. The numbers y^"^ are to the right of the vertical line. During 
each time step, the extremal site and its two neighbors are renewed by the triple of numbers to 
the right of the vertical line. 

figure. Each row is obtained by aligning the three y values of the row under the minimum 
X value of the previous row. 

B. P-completeness of standard dynamics 

P-completeness for a natural decision problem associated with standard Bak-Sneppen 
dynamics is proved by a reduction from the circuit value problem (CVP). 

The inputs of the STANDARD BAK SNEPPEN decision problem arc the initial site 
values, {xi(0)}, the renewal numbers, {y^"''}, a specified site k, time tmax and bound E. The 
problem is to determine whether Xk{tmax) < E. 

In order to effect the reduction from CVP we need to show how to build logical gates 
and wires that carry truth values. We will show how to build AND gates, OR gates and 
wires in the one-dimensional Bak-Sneppen model. This set of gates allows us to build 
monotone circuits. Since monotone CVP is P-complete only for non-planar circuits, our 
construction constitutes a P-completeness proof only after it is extended to > 1 Bak- 
Sneppen models. After we present the basic construction for o? = 1 we will discuss how to 
extend the construction to d> 1 and the prospects for proving P-completeness ior d— 1. 

In order to embed a circuit in standard dynamics, most of the sites of the lattice will have 
the value 1. A truth value is represented by a pair of neighboring sites i and such that 
within the pair, one site has the value one and the other has a value less than one. If the 
left member of the pair is a one the pair represents TRUE, otherwise it represents FALSE. 
A truth value is activated when its smaller member is the extremal site of the lattice. Once 
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the truth value is activated, it can be moved by appropriate choice of the triple of y's. The 
triple (1, 1, a) where a is a suitable chosen number, < a < 1 will move a truth value to the 
right while the triple (a, 1, 1) will move the truth value to the left. AND gates are realized by 
transporting two truth values next to one another. Suppose that initially, two truth values 
A and B are next to each other as shown in the Fig. |^. The numbers a, b, b', b" and c are 
ordered a < b < b' < b" < c and, at the time the gate is activated, c is less than all the other 
numbers in the lattice. The top panel in Fig. shows the situation where A = TRUE and 
B = TRUE. Three time steps later, C represented by the middle two sites of the lattice is 
TRUE. The top panel in Fig. shows the situation where A = TRUE and B = FALSE in 
which case C = FALSE as desired. The case A = FALSE and B = FALSE is the same as 
A = B = TRUE except the states in the part of the lattice representing the gate are shifted 
one lattice constant to the left so that C = FALSE as desired. The final case A = FALSE 
and B = TRUE is left as an exercise for the reader. 

OR gates are slightly more complicated than AND gates. Initially, the truth values are 
brought together until they are separated by one site as shown in Fig. The numbers are 
ordered < a < b < b' < b" < c < c' < c" < 1 and, at the time the gate is activated, c" is 
less than all the other numbers in the lattice. The top panel in Fig. ^ shows the situation 
where A = TRUE and B = TRUE. Six time steps later, C, represented by the second 
and third columns to the left of the vertical line, is TRUE. The situation A =FALSE and 
B =FALSE is seen to be correct by translational invariance. The case where A = FALSE 
and B = TRUE is shown in the bottom panel of Fig. ||. Verifying the case A = TRUE and 
B = FALSE is left as an exercise for the reader. 

We also need a gadget for fan-out. Figure |] shows how two copies, C and C can be 
obtained from A. The figure shows the case where A = TRUE. The case where A = FALSE 
is seen to be correct by translational invariance. 

Overall timing of the circuits is controlled by the numbers representing truth values and 
used in gates. Since, at each stage, the smallest number on the lattice is selected, truth 
values represented by small numbers are activated first. A truth value A that is active, 
meaning it is represented by a number a that is the minimum number on the lattice, can 
be put into storage for a pre-determined period by choosing y values 1) with b larger 
than all the values associated with logic variables that are active until A is active again. In 
this way all the required x and y values can be calculated with a fast parallel computation. 
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FIG. 7: Gadget for AND gate for standard dynamics. The inputs A and B are represented by 
the two pairs of sites above the horizontal hne. The output C is represented by the middle two 
sites along the bottom row to the left of the vertical line. The various numbers are ordered 
0<a<6<6'<6"<c<l. The top pictures shows the situation with both inputs are TRUE 
yielding an output of TRUE. The bottom picture shows the situation where the left input is TRUE 
and the right input is FALSE yielding FALSE. 

in advance knowing only the activation time of each logic gate. 

AND and OR gates, fan-out, timing and wires (transporting truth values) are sufficient 
for constructing arbitrary monotone Boolean circuits. Since planar monotone circuits are 
not P-complete, the preceding construction does not show that the d = 1 STANDARD 
BAK SNEPPEN problem is P-complete. It is straightforward to construct all the gadgets 
and arrange timings within d > 1 standard Bak-Sneppcn dynamics by simply padding the 
set of y values with I's in the directions transverse to the planes in which the action of 
the gadgets takes place. Since truth values can be transported in any direction, non-planar 
monotone circuits can be embedded in d > 1 Bak-Sneppen systems. The conclusion is that 
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FIG. 8: Gadget for OR gate for standard dynamics. The inputs A and B are represented by 
the two outer pairs of sites above the horizontal hne. The output C is represented by the second 
and third sites to the left of the vertical line. The various numbers are ordered < a < 6 < 
b' < b" < c < c' < c" < 1. The top pictures shows the situation with both inputs TRUE. The 
bottom picture shows the situation where the left input is FALSE and the right input is TRUE, 
also yielding TRUE. 

the d > 1 STANDARD BAK SNEPPEN problem is P-complete. Thus there can be no 
exponential speed-up of standard dynamics through parallelism (unless it should turn out 
that P = NC, a highly unhkely prospect). Although the current proof does not yield a 
P-completeness result for the planar STANDARD BAK SNEPPEN problem, we do not see 
a general reason preventing the construction of non-monotone circuits and we cannot rule 
out the possibility that the planar problem is also P-complete. 
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FIG. 9: Gadget for FAN OUT for standard dynamics. The input A is represented by the middle 
pair of sites above the horizontal line. The outputs C and C are represented by the two pairs of 
sites to left of the vertical line. The numbers are ordered a < c < c' . 

C. Power law speed-up of standard dynamics 

At least for > 1, the P-completeness result precludes the possibility of solving the 
STANDARD BAK SNEPPEN problem in polylog time. Nonetheless we can achieve a 
sublinear parallel time solution to the problem by invoking the hierarchical construction 
of avalanches from parallel dynamics. The history produced by standard dynamics can be 
viewed as a sequence of avalanches and these avalanches can be efficiently constructed from 



the list of y values using the parallel methods of Sec. [IV B| . The road block in achieving 
exponential speed-up for standard dynamics does not lie in the construction of avalanches, 
rather it is due to the presence of explicit initial values instead of initial probabilities defined 
by G{0) and made explicit only at the end of construction. The presence of explicit initial 
values means that the environment of each avalanche is specified and the sequence of gap 
values cannot be determined in advance of making the avalanches. 

The accelerated algorithm for the STANDARD BAK SNEPPEN problem proceeds as 
follows. The first gap value, G{0) = /(O) is just the minimum among the initial values. 
The algorithm then calls a parallel avalanche subroutine to produce a G'(0)-avalanche from 
the initial sequence of y values. This avalanche (the full avalanche, not just the backbone) 
renews the sites that it covers while the rest of the environment is unaffected. At the ending 
time, ti of the G'(0)-avalanche, the extremal site is determined and its values determines the 
next gap value, G(l) = f(ti). Generally, suppose that the gap increases at times t^, then 
G{k) = f{tk) and the parallel avalanche subroutine is called to produce a G(A;)-avalanche 
using the remaining list of y values. The origin of this avalanche is the extremal site s(tfc) 
and the avalanche renews all the sites it covers. As the gap approaches the critical value. 
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/c, the algorithm becomes increasingly efficient since larger and larger avalanches may be 
made in parallel. 

The parallel avalanche subroutine takes as inputs a gap value Gk and the remaining y 
values, {y^^^l n = {2d + + 1), . . . , (2d + l)tmax} values. Its output is the G^-avalanche 
that would have resulted if standard dynamics were used with the same sequence of y values. 
The procedure is similar to that described in Sec. |1V lj| and uses the same hierarchy of fn- 
avalanches but care must be taken that the y values are used in the correct order. Each 
block of 2d + 1 y values is an /o-avalanche. In parallel we must now group this sequence 
of /o-avalanches into /i- avalanches. The procedure for grouping avalanches is illustrated 
using Fig. ^. Each node on the bottom level of the tree represents blocks of 2d+l y values. 
In parallel we now independently start from each of these nodes and carry out standard 
dynamics until an /i-avalanche is complete. For example, starting from bottom level node 2 
it takes two time steps to build an /i-avalanche denoted by the 3 on the fi level of the tree 
but it only takes one time step starting from bottom level node 3 since this /o-avalanche is 
also an /i-avalanche. The time that standard dynamics must be run to obtain each of the 
/i-avalanches defining the nodes on the /i-level is just the maximum in-degree of this level. 
Thus far we have made /i-avalanches starting from each bottom level node, now we must 
determine which of these avalanches to keep to reproduce the result of standard dynamics 
acting on the y sequence. The correct /i-avalanche to associate with each node in the /i- 
level is the one that includes all of the children of the node. This choice insures that all y 
values are used. For example, the /i-avalanche labeled 3 includes steps 2 and 3 rather than 
3 alone. 

Having assembled the correct set of /i-avalanches we can now group these into /2- 
avalanches and so on. The new ingredient in going from the to /^-avalanches when 
n > 1 is that the generalized version of the standard dynamics is used as discussed in Sec. 



IV B| . As before, G^-avalanches can be constructed in a time that is polylog in the S, the 
duration of the avalanche. 

The speed-up obtained due to the parallel avalanche subroutine can be estimated by the 
scaling laws for the (S) and fc — {G(t)). Up to logarithmic factors, it takes constant parallel 
time to produce a G-avalanche. Combining Eqs. |1| and ^ we have 

dT 1 , 7 

(t/N) 7-1 (9) 



dt {S) 



G 
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where T is parallel time and t is sequential time. Integrating the differential equation yields, 

1 

T ti^^ SO that power law speed-up is obtained. 

The important difference between parallel and standard dynamics is that in parallel 
dynamics, sites that are not part of avalanches are specified probabilistically until the last 
stage of the computation while, for standard dynamics, all sites are explicitly determined 
on each time step. This suggests the idea that an efficiently parallelizeable dynamics might 
result from minimizing the number of explicitly determined values until the last stage of 
the computation. This idea lead to the conditional dynamics described and analyzed in the 
next two sections. 

D. Conditional Bak-Sneppen dynamics 

In standard dynamics with the usual choice of initial conditions, lattice sites are assigned 
random numbers uniformly distributed on [0, 1] and the lattice site with the smallest number 
is selected for replacement. In parallel dynamics only the backbone is fully specified until 
the end of the computation. In conditional dynamics this idea is carried as far as possible. 
Extremal space-time points are fully specified and conditions are placed on all other space- 
time points until the end of the calculation. Conditional dynamics is similar to the Run 
Time Statistics approach in Ref. |12|. 

Initially, each site is equally likely to be extremal. A site s(0) is selected at random 
from among all the sites of the lattice. Suppose there are N = L'^ sites on the lattice. The 
first extremal number, /(O) is chosen from the distribution of the minimum of numbers 
distributed uniformly on [0, 1]. Although we have not yet assigned numbers to any other 
sites we now know that these sites are conditioned to be greater than /(O). The next step in 
the standard Bak-Sneppen process is to renew s(0) and its neighbors with numbers chosen 
on the interval [0, 1]. Again, at this stage, we do not need to explicitly know any of the 
numbers on the lattice except to say that all the numbers on the lattice except s(0) and 
its neighbors are randomly chosen in the interval [/(O), 1] while s(0) and its neighbors are 
randomly chosen on the interval [0, 1]. We can characterize the condition on site i at time 
t by the number Cj(t) which indicates that the random number at site i is chosen from 
the uniform distribution on the interval [ci{t), 1]. We have that q(1) = for s(0) and all 
its neighbors while Ci(l) = /(O) for all other sites. The extremal site at t = 1 is a 
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weighted random choice among the N sites with site i weighted by 1 — Ci(l). The extremal 
value /(I) is chosen from the distribution of the minimum of numbers, where the i^^ 
number is chosen from the uniform distributions on [cj(l), 1]. The minimum numbers /(O) 
and /(I) now determine the conditions Cj(2) and so on. Specifically, s{t) and its neighbors 
are renewed at time t + 1 so that Ci{t + 1) = for s{t) and its neighbors. For all other sites, 
+ =max{/(t),Q(t)}. 

One way to implement the selection of the extremal site at time t given the conditions 
Ci{t) is to assign random numbers rj(t) on the interval [0, 1] to each space-time point. From 
these numbers, compute fj(t) = Ci(t) + ri(t){l — Ci(t)) = Ci(t) © ri{t). Note that ri(t) is 
uniformly distributed on [cj(t),l]. The extremal site s{t) is the site with the minimum 
among the f's and f{t) = rs{t){t)- Note that all the random numbers rj(t) can be assigned 
in advance. It is this realization of conditional Bak-Sneppen dynamics that is shown to be 
P-complete in Sec. [VI E| . 

After the conditional dynamics has run for a time tmax we are left with the conditions 
{ci{t)} together with a list of extremal sites {sit)\t < tmax} and extremal values {fit)\t < 
tmax}- From these data we can reconstruct a Bak-Sneppen history by working backwards. 
The approach is very similar to the construction of a full history from the backbone in 
the parallel dynamics. The final conditions, {ci{tmax)} can be used to obtain final values: 
Xiitmax) are randomly chosen from the uniform distributions on [ci{tmax)A]- For each i, 
these final values are then the correct values backwards in time through the latest time 
when site i is renewed. Suppose site i is last renewed at time t < tmax that is i = s(t — 1) 
or z is a neighbor of s{t — 1). If i = s{t — 1) then Xi{t — 1) = f{t — 1) and this value holds 
until the next earliest time that i is renewed. If site i is renewed at time t because it is the 
neighbor of the extremal site s{t — 1) then Xi{t — 1) is chosen from the uniform distribution 
on [max(/(t — l),Cj(t — 1)), 1] and this value holds until the next earliest time that i was 
renewed. Working backwards in this way a definite Bak-Sneppen history is constructed 
from the set of extremal values, extremal sites and conditions. It should be noted that the 
reconstruction of a definite history from the data generated by the conditional dynamics can 
be carried out in parallel in polylog time. 

During conditional dynamics, the values of the c's become spatially non-uniform however 
at times when the gap increases, almost all the c's are equalized. If G increases at time 
t — 1 then for all i, Ci{t) = G{t — 1) except, of course, at the extremal site s{t — 1) and its 
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neighbors where the c's are reset to zero. 



E. P-completeness of conditional dynamics 

P-completeness for conditional dynamics with > 1 is proved by a reduction from the 
monotone circuit value problem. To carry out the reduction we need a way to implement 
AND and OR gates and non-local fan-out. We will explain how the reduction works for 
the one-dimensional Bak-Sneppen model. Since non-local fan-out can be implemented, the 
P-completeness proof holds even for the one-dimensional problem. 

The inputs of the CONDITIONAL BAK SNEPPEN decision problem are the initial 
conditions, {ci(0)}, the random numbers, {ri(t)}, a specified site k, time tmax and bound E. 
The problem is to determine whether Ckitmax) < E. 

In the reduction from CVP, truth values reside at specified sites and are represented by 
the values of Ci{t). If Q(t) < et, the site has the value TRUE at time t and if Q(t) = 1 the 
site has the value FALSE. The number e is chosen to be small enough that the total time 
duration to evaluate the circuit, tmax still yields small values, thus we set etmax < 1/4. Sites 
with truth values are separated by background sites where rj(t) = Cj(t) = 1. If site i has an 
initial truth value this is represented by an initial value, Cj(0) = or 1. Sites that represent 
outputs of gates or fan-outs have Cj(0) = 1. 

In addition to sites representing truth values, there is a single site called $ with values 
r^{t) = et and c$(0) = 0. A site i carrying a truth value can be read at time t by setting 
rj(t) = 0. If site i is TRUE at time t it is selected since fj(t) < et whereas, if i is FALSE at 
time t then ri{t) = 1 and $ is selected at time t. 

Finally, for each AND gate and each fan-out, we need two neighboring auxiliary sites k 
and A with c^{0) = and ca(0) = 1. 

First, we show how to fan-out a truth value at site A at time t to a new site. A' at time 
t + 3 using the gadget shown in Fig. In Fig. |T0|, white squares have r = 1, black squares 



have r = and gray squares have small but non-zero r values. The gray squares in column 
$ have r$(t) = et and the gray square in column A at time t + 2 has rx{t + 2) = 3e/2. The 
value of Cyi(t) depends on the truth value of A. If A is TRUE then Cyi(t) < et but if A is 
FALSE then CA{t) = 1. At time t, site k is selected. First suppose A is TRUE. Then, since 
CA{t + 1) < et and ryi(t + 1) = we have that site A is selected at time t + 1. The selection 
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FIG. 10: Gadget for non-local fan-out of truth value at site A to site A'. 

of site A means that cx(t + 2) < et and therefore, fx(t + 2) < et + (3e/2)(l - et) < {t + 2)e 
so that site A is selected at time t + 2. The selection of A at time t + 2 insures that A' is set 
to TRUE for any time later than t + 2 since A' is refreshed at time t + 3 and ca' + 3) = 0. 
In Fig. [l^ we show A' read at time t + 3 but this is not necessary, it can also read at any 
time after t + 3. 

Now suppose A is FALSE. Then site $ is selected at time t + 1 and, consequently, 
f^(t + 2) = (t + l)e + (3e/2)(l - (t + l)e) > (t + 2)e. Thus site $ is selected at time t + 2 
and site A' fails to be refreshed so that CA'(t + 3) = 1 and A' is FALSE. 

The non-local fan-out just described can be easily extended to produce a non-local AND 



gate as shown in Fig. |TT]. As before, rA(t + 3) = 3e/2. The inputs A and B may be anywhere 
on the lattice, and the output AAB appears immediately to the right of the auxiliary sites 
K and A. Suppose first that both A and B are TRUE, then cx{t + 3) < et and fx{t + 3) < 
et + (3e/2)(l — et) < (t + 2)e so that A is selected at time t + 3 and AA B is refreshed at time 
t + 4. Thus the site AAB is correctly set to TRUE at time t + 4. On the other hand, suppose 
that A is FALSE and B is TRUE. Then fxit + 3) = (t + l)e + (3e/2)(l - (t + l)e) > (t + 3)e. 
Thus site AAB fails to be refreshed so that AAB is correctly set to FALSE. The other two 
possibilities for A and B are easily seen to work in the same way. 

The OR gate is local in space but non-local in time and is shown in Fig. |l^. Initially site 
AVB is set to FALSE however, at time t + 2, AVB is correctly set to A OR B since if either 
A or B are TRUE, then AVB is refreshed and can be selected at any later time to be the 
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FIG. 11: Gadget for non-local AND gate. Sites A and B may be anywhere in the lattice. At time 
t + 4, site AaB is set to A AND B. 
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FIG. 12: Gadget for OR gate. At time t + 3, site AvB is set to A OR B. 



input to a fan-out or AND gate. Figure shows A and B read at time t and t + 2 however 
any time separation between the reading A and B is permitted, the only constraint on the 
OR gate is that both A and B must be read before AVB is read. 

These gadgets are sufficient to reduce monotone CVP to the one-dimensional CONDI- 
TIONAL BAK SNEPPEN problem and show that the latter is P-complete. Note that 
the non-local character of the gadgets allows non-planar monotone CVP to be reduced to 
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the one-dimensional problem, in contrast to the situation for standard dynamics where the 
P-completeness holds only for d > 1. 
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VII. CONCLUSIONS 



Our main result is that Bak-Sneppen histories can be efficiently generated in paral- 
lel. Specifically, simulations and analytic arguments suggest that a history of length tmax 
can be simulated on a PRAM with polynomially many processors in average parallel time 
log'^*-^-' tmax with the actual asymptotic behavior close to log^ tmax • The exponential speed-up 
achieved by parallelization is the result of re-organizing the history into a nested hierarchy 
of independent avalanches. The construction of a single /-avalanche can be carried out in 
average parallel time log^*^^^(/c — /) which is exponentially less than the expected duration 
of the avalanche {S) j ^ {fc — f)~'^ ■ 

How is it possible to create long range correlations in space and time very quickly in 
parallel? First of all, it must be emphasized that the ground rules for parallel computation 
with a PRAM allow for non-local transmission of information in a single time step. To 
see how space-time correlations can be set up in a parallel time that is polylog in the 
correlation length or time, consider one step in the hierarchical construction of avalanches. 
For example, suppose two /„_i-avalanches are concatenated to yield a single /„-avalanche. 
The correlation length and time are thereby increased by constant factors in a fixed number 
of parallel steps. Non-local transmission of information is needed to align the origin of the 
second /„_i-avalanche with the final extremal site of the first /„_i-avalanche and it is this 
alignment of the two /„_i-avalanches that increases the correlation length and time. The 
result is that spatio-temporal correlations grow exponentially in the number of parallel steps. 

In the parallel construction, an mdependent collection of /„_i-avalanches are concatenated 
to form an /„-avalanche. Nonetheless, Bak-Sneppen avalanches have a temporal structure 



and exhibit aging |T^. Aging is consistent with the independence of subavalanches for the 
following reason. The last extremal value in the construction of an /^-avalanche must exceed 
fn and this extremal value is more likely to come from the last /„_i-avalanche used in the 
construction. Thus, extremal values and other properties of avalanches display aging. 

Another example where critical correlations are set up much faster in parallel than might 
be expected are cluster Monte Carlo algorithms |T5[ for critical spin systems. Each cluster 
sweep can be accomplished in polylog parallel time on a PRAM with polynomially many 
processors using a parallel graph connectivity algorithm. The number of sweeps required 
to reach equilibrium scales as a small power (typically much less than unity) of the system 
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size. Thus long range critical correlations are set up in a parallel time that is much less than 
the correlation length. This is quite different than physically realistic local dynamics where 
critical correlations require a time that is at least quadratic in the correlation length. 

The intuition that Bak-Sneppen histories are generated by an inherently sequentially 
process is not entirely wrong. In fact, standard Bak-Sneppen dynamics was shown here to 
be associated with a P-complete decision problem for o? > 1. Standard dynamics takes as its 
inputs the initial values on the lattice and a sequence of {2d + l)-tuples of numbers that are 
used to renew the extremal site and its neighbors at each time step. The P-completeness 
result implies that it is almost certainly not possible to parallelize standard dynamics so 
that the resulting history is generated in polylog time. The key problem that prevents full 
parallelization of standard dynamics is that all site values arc explicitly defined at every 
step in the construction. In contrast, parallel dynamics does not assign explicit values 
except in the avalanche backbones until the end of the calculation. For the Bak-Sneppen 
model, predicting the future starting from explicit initial conditions is harder than sampling 
a typical history. 

Conditional dynamics is a third method for generating Bak-Sneppen histories that has the 
smallest set of space-time sites explicitly defined until the end of the calculation. Nonethe- 
less, conditional dynamics is also associated with a P-complete decision problem and so 
cannot be efficiently parallelized. Parallel dynamics shares features in common with both 
conditional and standard dynamics and is intermediate in the degree to which it avoids 
explicit specification of site values until the end of the computation. 

The results presented here refiect characteristics of the Bak-Sneppen model that are 
independent of the PRAM model of computation in which they were presented. The pri- 
mary result of the paper can equivalently be stated in terms of families of Boolean circuits. 
Specifically, we can envision a hard-wired device composed of logical gates and random bit 
generators. The gates are wired in a feedforward direction so that, in generating a history, 
each gate is used only once. When the circuit is activated it produces a statistically correct 
Bak-Sneppen history. For each system size L and time tmax we need a different circuit. Be- 
cause of the equivalence of PRAMs and circuit families, our result can be stated in terms of 
the logical depth of the circuit. Logical depth is the longest sequence of gates between any 
of the random bit generators and the output representing the history. The existence of an 
efficient parallel algorithm implies that there is a uniform family of the circuits for generat- 
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ing Bak-Sneppen histories whose depth scales polylogarithmically in tmax, independent of L. 
The number of gates in the circuit is bounded by a power of t^^^^L. The actual running time 
for any real circuit generating Bak-Sneppen histories would, of course, be polynomial and 
not polylogarithmic in tmax because of the need for wires connecting distant gates required 
for establishing the long range spatio-temporal correlations. 

A variety of non-equilibrium models in statistical physics have been employed to study 
the spontaneous emergence of complexity. Besides the Bak-Sneppen model, other examples 
include diffusion limited aggregation ||l6l (DLA), invasion percolation and sandpiles. Besides 
their application to specific physical phenomena, these models have a broad appeal because 
they are governed by simple microscopic rules and yet they display self-organized criticality. 
It has been argued that they perhaps shed light on far more complex phenomena found, 
for example, in biology or economics. Surely, however, biological and economic systems 
generate histories that have polynomial rather than polylogarithmic logical depth. Perhaps 
one criterion for a model of the spontaneous emergence of complexity is that it should 
generate histories that require more than polylogarithmic depth to simulate. Neither the 
Bak-Sneppen model nor invasion percolation |^ satisfy this requirement. On the other 
hand, both diffusion limited aggregation and the Bak, Tang and Wiesenfeld (BTW) sandpile 
model are related to P-complete problems |^ which strongly suggests that the clusters 
or avalanches associated with these models cannot be simulated in polylog depth. However, 
neither model generates histories much deeper than a power of the system size. In the case 
of DLA, the cluster size cannot exceed the system size so that the length of the history is 
bounded by tmax < where D is the fractal dimension of the cluster. Sandpile models are 
well-defined for arbitrarily long times. However, the BTW model has an Abelian property, 
which allows avalanches to be combined in any order. Thus, it is possible to rearrange long 
histories so they are generated in polylog parallel time. From a given initial condition, it 
requires C(L'^+^) steps to relax the system 0. A sequence of avalanches due to randomly 
dropping sand on single sites can be organized into a binary tree using the same idea as 
adding numbers in parallel. The result is that a history of length tmax can be simulated in 
parallel in time log tmax)- By contrast, we find that Bak-Sneppen histories can be 

simulated in parallel time log'^*^^-* tmax independent of system size. 
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